AD-A145  818  COMPUTER  MOLECULAR  DVNAMICS  OF  2-D  WATER  (ICE) 

STRUCTURES  UNDER  SHOCK(U)  LAWRENCE  LIVERMORE  NATIONAL 
LAB  CA  P  HARRIS  ET  AL.  AUG  84  UCID-20137 
UNCLASSIFIED  N08014-84-F-0094  F/G  8/12 


AD- A 145  818 


IIC  ID-20137 


COMPUTER  MOLECULAR  DYNAMICS  OF  2-D  WATER  (ICE) 
STRUCTURES  UNDER  SHOCK 
-  A  STATUS  REPORT  - 


Paul  Harris 
Arnold  M.  Karo 


VV  ot?isee4 


August  1984 


This  is  an  informal  report  intended  primarily  for  internal  or  limited  external 
distribution.  The  opinions  and  conclusions  stated  are  those  of  the  author  and 
may  or  may  not  be  those  of  the  Laboratory. 

Work  performed  under  the  auspices  of  the  U.S.  Department  of  Energy  by  the 
Lawrence  Livermore  National  Laboratory  under  Contract  W-7405-Eng-48. 


r<& 


*5 


84  08  28  010 


DISCLAIM  KR 


This  document  was  prepared  as  an  account  of  mirk  sponsored  by  an  agency  of  the  I  nited  States  (Government. 
Neither  the  llnited  States  (Government  nor  the  I  Diversity  of  (  alifornia  nor  any  of  their  employees,  makes 
any  warranty,  express  or  implied,  or  assumes  any  legal  liability  or  responsibility  for  the  accuracy,  com¬ 
pleteness,  or  usefulness  of  any  information,  apparatus,  product,  or  process  disclosed,  or  represents  that  its 
use  would  not  infringe  privately  owned  rights.  Reference  herein  to  any  specific  commercial  products,  process, 
or  service  by  trade  name,  trademark,  manufacturer,  or  otherwise,  docs  not  necessarily  constitute  or  imply  its 
endorsement,  recommendation,  or  favoring  bv  the  t  nited  States  (Government  or  the  I  Diversity  of  (  alifornia. 
The  views  and  opinions  of  authors  expressed  herein  do  not  necessarily  slate  or  reflect  those  of  the  I  nited 
States  Government  thereof,  and  shall  not  be  used  for  advertising  or  product  endorsement  purposes. 


Printed  in  the  I'm  led  Stales  of  America 
Available  from 

National  technical  Information  Service 
U.S  Department  of  Commerce 
■>28^  Port  Koval  Koad 
Springfield,  VA  22 IM 

Trite:  Printed  i.  npv  S  ;  Mu'mliche  S4 -*>ll 


Page  Range 

Domestic 

Price 

001-025 

$  7.00 

026-050 

8.50 

051-075 

10.00 

076-100 

11.50 

101-125 

13.00 

126-150 

14.50 

151-175 

16.00 

176-200 

17.50 

201-225 

19.00 

226-250 

20.50 

251-275 

22.00 

276-300 

23.50 

301-325 

25.00 

Page  Range 

Domestic 

Price 

326-350 

$  26.50 

351-375 

28.00 

376-400 

29.50 

401-426 

31.00 

427-450 

32.50 

451-475 

34.00 

476-500 

35.50 

501-525 

37.00 

526-550 

38.50 

551-575 

40.00 

576-600 

41.50 

bOl-up' 

’Add  1.50  for  each  additional  25  page  increment,  or  portion 
thereof  from  601  pages  up. 


COMPUTER  MOLECULAR  DYNAMICS  OF  2-D  WATER  (ICE) 
STRUCTURES  UNDER  SHOCK 
-  A  STATUS  REPORT  - 

by 


Paul  Harris 

US  Army  Armament,  Munitions  and  Chemical  Command 
Dover,  New  Jersey  07801 


Arnold  M.  Karo 

Lawrence  Livermore  National  Laboratory 
Livermore,  California  94550 


ABSTRACT 


Computer  molecular  dynamics  has  been  used  to  study  the  shock-front 
transition  region  for  a  series  of  two-dimensional  hydrogen-bonded 
structures  composed  of  water-vapor  molecules.  Two  different  structures 
were  examined  and  different  intermolecular  potentials  were  considered. 
Sample  sizes  and  running  times  were  chosen  to  correspond  to  the  predicted 
shock-front  rise-time  in  real  water.  In  this  way  the  effect  of  different 
potentials  and  initial  structures  on  the  equilibration  associated  with 
the  shock-front  transition  in  water  could  be  investigated.  In  addition 
to  studying  the  development  and  incipient  relaxation  of 
shock-polarization  states,  we  have  also  considered  the  propagation  and 
possible  incipient  relaxation  of  structural  phase  transitions  occuring  in 
two-dimensional  structures.  The  results  of  these  molecular  dynamics 
calculations  are  compared  with  the  experimental  shock  behavior  of  real 
water;  in  particular,  comparisons  are  made  with  respect  to  intermolecular 
hydrogen-bond  breaking,  dissociation,  and  dissociation-related  electrical 
conductivity  data.  », 


0R52M/dsm 


INTRODUCTION 


The  problem  of  providing  a  detailed  (although  averaged)  description 

of  molecular  behavior  within  a  shock-front  transition  region  is  of  both 

academic  and  applied  interest.  The  academic  interest  is  that,  except  for 

1-3 

the  case  of  a  condensed  monatomic  noble  gas  ,  the  problem  has  never 
been  successfully  treated  for  realistic  condensed  molecular  media.  The 
applied  interest  is  mostly  connected  to  understanding  the  possibility  of 
enhanced  cnemistry  occurring  within  the  high-strain-rate  shock-front 
region. 

Water  was  chosen  for  this  study  because  water  has  been  extensively 
studied  both  macroscopically  and  microscopically;  because 
shock-polarization4’^  and  shock-front  optical  reflectivity^  data 
exist  for  water;  and  because,  to  a  first  approximation,  each  molecule 
could,  if  desired  for  purposes  of  simplicity,  be  considered  as  a  rigid^ 
structure  in  the  presence  of  an  intermolecular  potential.  Thus,  water 
should  represent  a  benchmark  test  of  the  ability  of  computer  molecular 
dynamics  (CMD)  to  predict  the  known  shock -wave  properties  of  a  realistic 
and  interesting  material. 

The  work  to  be  described  in  this  report  represents  the  first  steps 
along  the  road  to  a  full  and  accurate  CMD  study  of  shocked  water.  Two 
different  initial  two-dimensional  (2D)  structures  were  run.  One  of  these 
structures  was  also  run  with  important  differences  in  the  intermolecular 
potentials.  The  problems  were  run  only  in  2D  configurations  because  the 
most  important  aspect  of  the  physics  (chemistry)  of  water  is  that 
associated  with  the  intermolecular  hydrogen  bonding,  which  can  be 
adequately  represented  in  20,  and  to  save  time  (as  well  as  expense) 
during  the  actual  computations.  In  our  calculations  we  have  not  found  it 
necessary  to  treat  the  water  molecular  as  rigid  and  have  therefore 
allowed  for  all  2D  molecular  degrees  of  freedom. 

Among  the  shock-related  phenomena  studied  were  the  propagation  and 
possible  commencement  of  relaxation  of  structural  phase  transitions,  the 
development  and  possible  commencement  of  relaxation  of  shock-polarization 
states,  and  the  importance  of  the  details  of  the  potentials  to  the 
general  questions  of  realistic  shock  wave  CMD. 
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For  these  studies  the  largest  impacted  sample  considered  was  composed 
of  4  x  20  water  molecules  (20  molecules  in  the  direction  of  shock 
propagation),  and  the  longest  running  time  (measured  from  impact)  was 
0.65  picoseconds.  These  sample  length  and  run  times  correspond 
approximately  to  the  predicted4  shock-front  rise  time  in  real  water  of 
approximately  one  picosecond.  Thus,  at  least  in  principle,  the  effect  of 
potentials  and  initial  structures  upon  the  equilibration  associated  with 
the  shock-front  transition  in  water  could  begin  to  be  investigated  with 
the  CMD  techniques  presented  in  this  report. 

This  report  presents  results  only  for  impact  velocities  of  0.75  A 
per  tick  (or  7.5  x  10  cm/sec)  and  0.4  A  per  tick.  Such  impact 
velocities  represent  (for  real  3-dimensional  water  under  1-dimensional 
strain  with  symmetric  impact)  pressures  of  approximately  300  kbar  (or 
30  GPa)  and  100  kbar  respectively. 

a 

The  electrical  conductivity  of  water  is  known  from  experiments  to 
saturate  at  shock  pressures  above  300  kbar.  It  is  thus  possible  to 
characterize  real  water  above  300  kbar  as  being  entirely  ionized. 

Ionized  water  has  been  described  as  consisting  primarily  of  protonated 
water  complexes10  (e.g.,  Hg0^)  and  hydroxide  ion  radicals.  In 
so  far  as  such  protonated  complexes  or  hydroxide  radicals  were  not 
observed  (i.e.,  no  recognizable  "dissociation”)  the  simple  CMD  model  we 
have  employed  cle?rly  requires  additional  refinement  in  order  to  be  in 
agreement  with  this  picture.  Thus,  at  this  moment,  although  it  is  not 
completely  clear  what  the  difficulty  is,  an  improved  model  must  include 
the  pressure  or  density  dependence  of  the  potential  barrier  to  proton 
transfer.  It  appears  that  although  the  shock  contains  sufficient  energy 
to  overcome  the  original  tetrahedral  orientation  and  average  bonding 
configuration  of  the  water  molecular  lattice  (as  discussed  in  the 
following  sections),  because  of  the  absence  of  energetically  low-lying 
pathways  that  energy  is  not  being  transferred  during  the  time  of  our 
simulation  studies  into  the  water  molecule  vibrational  modes  leading  to 
proton  transfer  and  the 


+  One  "tick" 
CMD  (10'14 


is  10  femtoseconds  or  10“14  sec  -  a  convenient  time  in 
sec  is  approximately  one  vibrational  period). 
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formation  of  these  complex  ionic  species.  (In  our  view  the  lack  of 
observed  complex  formation  at  the  0.75  A/tick  impact  velocity  indicates 
clearly  a  direction  for  further  study.) 

With  the  exception  of  the  intermolecular  hydrogen-hydrogen 
interactions,  all  of  the  interatomic  potentials  employed  were  of  the 
Morse  type  (with  a  cut-off  distance  and  damping  factor).  Some  of  the 
intermolecular  hydrogen-hydrogen  potentials  were  chosen  to  be 
Coulomb-like  (based  upon  the  molecular  orbital  hydrogen-bond  studies  of 
Morokuma  and  Pedersen^). 

The  above-mentioned  choice  of  intermolecular  potentials  describes 
firstly  a  relatively  fixed  (i.e.,  a  steep-sided  moderately-shal low 
potential  well)  OH  bond  distance  with  a  minimum 'ul  .76  A.  Secondly, 
the  choice  describes  the  separation  of  water  molecules  which  are  hydrogen 
bonded  to  a  common  oxygen  atom.  That  second  point  leads  to  a  very 
shallow  and  gradually  sloping  00  potential  well  yielding  a 
nearest-neighbor  00  distance  of  ^2 .76  A. 

With  the  above  OH  and  00  intermolecular  potentials,  together  with  a 

reasonable  description  of  the  water-vapor  molecule,  it  is  not  necessary 

to  require  a  potential  well  for  the  HH  intermolecular  interaction  because 

the  basic  water  structure  (in  the  sense  of  an  ice)  is  approximately 

fixed.  Accordingly  either  a  Coulomb  or  a  Morse  type  of  interaction  can 

be  chosen  for  the  intermolecular  HH  interaction.  This  choice  allows  for 

a  parameter  (charge  in  the  case  of  the  Coulomb  interaction)  which  when 

varied  will  have  an  influence  upon  the  ease  of  rotation  of  one  water 

molecule  in  the  field  of  another.  The  rotational-influence  point-of-view 

is  also  in  keeping  with  our  desire  to  ultimately  correlate  the  CMO  work 

4  5 

with  existing  shock-polarization  data  ’  (shock  polarization  is 
intimately  connected  to  the  shock-induced  average  rotation  of  the  water 
molecule  electric  dipole  moment). 

The  CMD  calculations  were  carried  out  with  initially  quiescent  (T  =  0 
degree  Kelvin)  lattices.  The  basic  reason  for  the  T  =  0  K  choice  is  that 
for  a  Hugoniot  pressure  of  ^300  kbar  the  shock-induced  temperature 
change  (^2000  K)  is  so  large  compared  to  room  temperature  that  it  is 
not  likely  that  an  initial  temperature  of  ^300  K  could  have  any  effect 
upon  the  qualitative  results  being  investigated  here.  With  that 
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understanding,  any  initial  thermal izat ion  of  the  lattice  would  have 
represented  an  unnecessary  expenditure  of  time  and  money  for  these 
first-step  calculations. 

All  of  the  CMD  problems  were  run  with  periodic  boundary  conditions  in 
the  direction  transverse  to  the  direction  of  shock  propagation.  This 
means  that  if  an  atom  is  at  position  (x,  y),  then  there  are  also 
identically  behaving  atoms  at  positions  (x,  y  +  nL)  for  all  integer  n 
where  L  is  the  transverse  sample  dimension  "size".  To  insure  basic 
structural  stability,  each  sample  (with  related  potentials)  was  "run" 
without  periodic  boundary  conditions  (in  the  absence  of  impact).  No 
phase  transitions  were  noted  in  those  stability  tests,  and  where  internal 
heating  occurred  due  to  small  initial  intrinsic  strain  it  was  always 
correctable  by  reducing  the  error  in  the  initial  atomic  equilibrium 
positions. 

II.  SOME  PRELIMINARY  NUMBERS 


Before  proceeding  with  the  detailed  presentation  of  the  actual 
potentials  and  structures  employed  in  these  CMD  calculations  it  is 
worthwhile  to  discuss  some  order-of -magnitude  numbers  relating  to  real 
water  and  its  behavior  under  shock. 

A.  An  average  value  for  the  water-vapor  dissociation  energy 

(removal  of  one  hydrogen)  is  approximately  120  kcal/mole  or 
-IP 

8.33  x  10  ergs/molecule. 

B.  The  symmetrical  bending  mode  for  a  water-vapor  molecule  has  a 

13  - 1 

characteristic  frequency  f  of  approximately  1600  cm  . 

2 

For  a  spring  constant  K  given  by  K  =  mH  (2irf)  with  m^ 

being  a  hydrogen  mass,  there  is  an  associated  bending  energy  of 
2 

E  =  0.5  K  (AX)  ,  where  AX  is  the  hydrogen-hydrogen 

vibrational  amplitude.  It  can  be  seen  that  a  (AX)  value  of 

-13 

0.2  A  corresponds  to  approximately  3  x  10  ergs. 

C.  The  hydrogen-bond  (intermolecular )  energy  responsible  for  the 

14 

basic  structure  of  water  and  ice  is  approximately 

-13 

5  kcal/mole  or  3.4  x  10  ergs/molecule . 


At  10  kbar  each  water  molecule  acquires  a  particle  (flow) 
velocity  (4.26  x  104  cm/sec)  kinetic  enerqy  of 
0.5  m^Up  or  2.6  x  ergs.  mH^  is  the 
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water-vapor  molecular  mass  and  u  the  particle  velocity. 

P  O 

There  is  also  a  PdV  contribution  (for  a  20%  strain  )  of 

approximately  the  same  energy  per  molecule.  The  water 

temperature  is  raised®  approximately  34  K  to  327  K  so  that 

after  equilibration  each  degree  of  freedom  is  characterized  by  a 

-14 

thermal  energy  of  0.5  kT  or  2.4  x  10  erg. 

E.  At  100  kbar  the  flow  and  PdV  contributions  per  molecule  are  each 

-13 

approximately  5.7  x  10  ergs  while  each  degree  of  freedom  is 

-14 

characterized  by  843  K  or  5.8  x  10  erg. 

F.  At  300  kbar  the  flow  and  PdV  contributions  per  molecule  are  each 

-12 

approximately  2.2  x  10  ergs  while  each  deqree  of  freedom  is 

-13 

characterized  by  2310  K  or  1.5  x  10  erq. 

From  A  through  F  above,  the  following  is  to  be  expected  with 
increasing  shock  pressure.  At  10  kbar  neither  intermolecular  OH  bond 
breakinq,  intramolecular  HH  separation  changes,  or  dissociation  should 
occur.  At  100  kbar  intermolecular  OH  bond  breaking  and  changes  in  HH 
intermolecular  separation  should  occur,  but  there  should  not  be  any 
water-vapor  molecule  dissociation.  The  100  kbar  expectations  also  hold 
at  300  kbar. 

What  is  seen  in  the  CMD  runs  to  be  discussed  below  is  intermolecular 
OH  bond  breaking  and  intramolecular  HH  separation  distance  changes  at  100 
kbar  and  above.  In  accordance  with  the  binding  energies  associated  with 
the  potentials  we  are  using,  no  water-vapor  molecule  dissociation  is 
seen.  Thus,  with  the  exception  of  the  previously-discussed  ionization 
via  proton  transfer  interpretation  of  the  experimental  electrical 
conductivity  saturation  near  300  kbar,  the  CMD  runs  are  behaving  as 
desired. 

Thus  the  various  intermolecular  and  intramolecular  energies  are  a 
benchmark  not  only  for  the  potentials  to  be  described  in  the  next 
section,  but  also  a  benchmark  for  the  consequent  behavior  of  the  shocked 
CMO  water  structures. 


t 


» 
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III.  THE  POTENTIALS  AND  STRUCTURES 

With  the  exception  of  the  HH  intermolecular  potentials,  all  of  the 
potentials  employed  were  of  the  Morse  type: 

f  -A,(R-R  )  I  2 

V (R)  =  0e  I  e  1  ~1J  (1) 

where  L 

Rg  =  The  equilibrium  interatomic  separation  (where  V ( R )  achieves  a 
minimum) . 

-Dg  =  The  energy  of  the  potential  minimum  relative  to  V(<»). 

A-j  =  The  small  displacement  (relative  to  Rg)  force  constant 
parameter  with  the  small  displacement  force  given  by 
-2DeA^(R-Re). 

The  HH  intermolecular  potentials  were  either  of  the  Morse  type, 

Eq.  (1),  or  the  Coulomb  type,  given  by 
n2 

V(R)  =  a  (2) 

where 

a  =  The  attraction-repulsion  parameter.  a=  +1  for  repulsion 
or  -1  for  attraction. 

Q  =  The  interaction  strength. 

Table  I  gives  a  summary  of  the  three  CMD  runs  to  be  presented  in  this 
report  in  terms  of  intermolecular  potentials,  impact  velocities,  and  the 
starting  lattice  structures,  figures  1  and  2  illustrate  the  two  starting 
structures  (henceforth  called  triangular  and  square,  respectively)  with 
the  related  initial  distances  and  angles. 

Aside  from  the  difference  in  unit-cell  configurations  between  the  two 
starting  structures  (and  a  corresponding  compressibility  difference),  the 
lattices  are  characterized  by  different  responses  to  local  rotation 
because  of  the  different  hydrogen  arrangements  and  the  associated 
intermolecular  HH  potentials.  Such  rotational  response  differences  are 
important  to  the  question  of  shock  polarization. 
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With  the  exception  of  the  HH  intermolecular  Coulomb  potential,  the 
numerical  details  of  the  potentials  are  given  in  Tables  II  and  III.  The 
Coulomb  intermolecular  HH  potentials  used  in  runs  2  and  3  are  illustrated 
in  Fig.  3.  Other  than  for  the  cut-off  distances,  the  potential 
functional  forms  were  the  same  with  a  =  1,  Q  =  0.139  e  where  e  is  the 
electron  charge.  The  Q  =  0.139  e  value  arises  from  requiring  that,  at 
equilibrium  separation, 

n2 

-D  (00)  -  D  I  OH )  +  =  5  kcal/mole.  (3) 

e  e  Ke 

The  same  value  of  Q  is  conceptually  justified  for  use  in  both  runs  2  and 
3  by  the  fact  that  the  value  of  Q,  although  sensitive  to  Dg( OH )  and 
De(00),  is  rather  insensitive  to  the  HOH  spatial  orientation  relative 
to  the  0H0H-,  plane  as  shown  by  Morokuma  and  Pederson11.  Indeed  that 
insensitivity  can  be  viewed  as  a  justification  for  modeling  water  in  two 
dimensions. 

All  of  the  potentials  used  in  these  CMD  runs  contained  a  "damping" 
term  which  smoothed  the  potential  in  its  approach  to  a  constant  value  at 
the  cut-off  distance  Rc.  This  was  done  in  order  to  prevent  a 
discontinuity  in  the  forces  at  Rc.  For  example,  run  1  utilized  a 
damping  function  which  led  to  a  smoothed  force,  Fp,  given  by 

F0(R)  -  F(R)  -  F(RC)  exp  -16  (^)  ]  (4) 


Rc  -  R 
Rc  -  Re 


where  F(R)  is  the  force  derived  from  the  undamped  potential. 

The  magnitudes  of  the  potentials  for  the  different  intermolecular 
interactions  (all  runs)  at  equilibrium  (net  equilibrium  in  the  Coulomb 
interaction  cases)  are  such  that 


|V(HH)|  «  IV ( OH )l  ~  10| V (00)|  (5) 

While  potential  well  depth  is  important  for  "chemistry"  (e.g.,  bond 
breaking),  the  dynamics  are  controlled  by  the  curvature  of  the  potentials 
near  equilibrium.  For  example,  for  a  0.1  A  displacement  from 
equilibrium  in  the  run  3  Coulomb-interaction  case,  the 
intermolecular-force  magnitude  ratios  are  given  by  (from  Table  III  and 
Q  =  0.139) 
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(6) 


The  inter-  and  intramolecular  forces  must,  of  course,  be  related  to 
known  water  inter-  and  intramolecular  vibrational  frequencies.  From 

p 

F  =  -20  Af(R-R  )  it  is  seen  that  near  equilibrium  a  Morse 

C  I  C  A 

potential  is  equivalent  to  a  spring  constant  K  =  2DeA^  with  a 
corresponding  frequency  f  (in  inverse  cm  optical  units)  of 


2Defll 


(7) 


where  c  is  the  speed  of  light  in  vacuum.  Using  the  reduced  mass  y  for 
OH  in  Eq.  (7)  gives,  from  Table  II,  for  the  intramolecular  OH  interact 


fintra<0H>  ’  3'9  *  '°3 


(8a) 


From  Table  III  and  with  the  reduced  mass  corresponding  to  the  interaction 
of  two  water  molecules. 


finter(°H)  =  1 *9  x  1q3  cnfl  *  run  1  (8b) 

finter(°H)  =  °*4  x  1q3  cm-1,  runs  2  and  3.  (8c) 


The  Eq.  (8a)  re_,ult  corresponds  fairly  close  to  reality  in  that  it  is 
a  frequency  approximately  equal  to  that  associated  with  the  water-vapor 
molecule  symmetrical  zero-order  stretching  mode  (3825.3  cm”  )  . 

Eqs.  (8b)  and  (8c)  represent  a  hard  shallow  potential  (from  Table  III, 
0.57  eV  for  Eq.  (8b)  and  0.52  eV  for  the  Eq.  (8c)  cases).  There  is  no 
single  piece  of  experimental  evidence  with  which  to  compare  the  Eqs.  (8b) 
and  (8c)  results  in  that,  for  example,  the  observed^  (Raman  spectra) 
intermolecular  hydrogen  bond  stretching  mode  frequency  of  approximately 
200  cm'^  corresponds  to  the  totality  of  potentials  which  via  their 
interaction  results  i  the  "hydrogen  bond"  within  the  total  water 
structure. 

The  physical  effect  of  the  hard  (Eqs.  (8b)  and  (8c))  intermolecular 
OH  potentials  used  in  these  CMD  runs  is  that  OH  bond  breaking  will  occur 
for  a  smaller  OH  bond  length  change  than,  say,  for  a  potential  whose  bond 


frequency  is  characterized  by  200  cm" ^ .  The  energy,  however,  for  the 
proton  to  be  raised  out  of  its  well  (bond  breaking)  would  be  the  same 
('vO.S  eV)  for  either  the  Eq.  (8b)  and  Eq.  (8c)  frequencies  or 
200  cnf^ .  However,  as  indicated  earlier,  under  the  pressures  and 
densities  characteristic  of  the  shock  front  and  the  region  behind  the 
front,  the  potential  barrier  to  proton  transfer  could  be  considerably 
reduced  or  eliminated.  This  effect  has  not  been  considered  in  the 
present  calculations  but  certainly  should  be  included  at  a  later  stage. 

Because,  as  seen  from  the  Eq.  (6)  force  ratios,  the  hardest 

intermolecular  potential  is  associated  with  the  OH  bonding  potential  it 

is  possible  to  interpret  the  above  paragraph  in  a  shock-compression  PdV 

language.  There  are  two  intermolecular  OH  bonds  per  molecule  and 

22  T 

approximately  3.5  x  10  water  vapor  molecules  per  cm  .  For  our 
present  pressure-independent  intermolecular  potentials  the  0.52  eV  OH 
bonding  wells  are  equivalent  to  a  potential-energy  pressure  of 
approximately  58  kbar.  Assuming  that  the  shock  pressure  is 
equipartitioned  between  particle  velocity  and  strain  one  would  thus 
expect  significant  intermolecular  bond  breaking  in  these  CMD  runs  in  the 
vicinity  of  120  kbar  (as  indeed  observed  in  run  3).  If  the  OH  bond 
strength  was  softer  and  more  competitive  with,  say,  the  00  intermolecular 
potential,  then  the  above  shock-pressure  interpretation  would  be  changed 
(to  higher  pressures)  because  an  increased  fraction  of  the  total  strain 
energy  would  then  be  going  into  the  00  bonds+. 

The  2.05  A  cut-off  distance  for  the  intermolecular  HH  potentials 
with  the  triangular  structure  is  smaller  than  any  initial  inter-row  HH 
separations.  Thus  for  the  triangular  structures  (runs  1  and  2)  molecular 
rotation  transfer  across  molecular  rows  is  expected  to  be  inhibited  -  as 
indeed  observed  in  the  CMD  results.  On  the  other  hand,  the  distance 
between  atoms  H^  and  Hb  in  the  Fig.  2  square  lattice  is  less  than  the 

+  Because  an  00  bond  is  col  inear  with  a  very  hard  intramolecular  OH  bond 
plus  a  hard  intermolecular  OH  bond,  the  00  and  OH  intermolecular  bonds 
are  subject  to  the  same  change  in  length  for  a  compression  along  their 
common  bond  direction. 
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run  3  cut-off  distance  of  2.3  A.  Thus  for  the  square  lattice  one  would 
expect  to  see  a  more  uniform  (row  to  row)  molecular  rotation  distribution 
pattern  than  in  the  triangular  lattice  case  -  again  this  is  borne  out  in 
the  CMD  runs. 


IV.  THE  CMD  RESULTS 

Molecular  dynamics  has  now  evolved  to  a  sophisticated  level  because 
of  the  dramatic  revolution  in  large-scale  computer  technology.  Simply 
stated,  molecular  dynamics  involves  the  numerical  solution  by  computer  of 
Newton's  equations  of  motion  for  all  the  atoms  comprising  the  active 
region  of  the  assembly.  As  a  result  the  coordinates  and  velocities  of 
the  particles  are  obtained  as  functions  of  time. 

The  force  acting  on  any  one  particle  in  the  ensemble  is  the  resultant 
of  all  interactions  with  other  atoms  in  its  neighborhood.  In  principle 
it  should  be  obtained  as  the  derivative  of  an  effective  many-body 
potential.  However,  any  calculations  which  would  in  fact  attempt  to 
include  most  or  all  of  a  complete  potential  energy  hypersurface  would  be 
far  beyond  present  considerations.  Fortunately,  our  experience  to  date 
indicates  that  much  can  be  learned  within  the  framework  of  the  simple 
pair-potential  approximation. 

The  set  of  cojpled  non-linear,  second-order  differential  equations 
that  represent  the  forces  are  reduced  to  a  double  set  of  first-order 
ordinary  differential  equations  by  standard  procedures.  The  initial 
positions  and  velocities  of  the  particles  represent  the  initial 
conditions  on  these  equations.  For  the  numerical  solution  of  this  set  of 
equations  we  have  employed  the  variable-step,  variable-order,  implicit 
Adams'  method  with  functional  (or  fixed  point)  iteration;  a 
generalization  of  the  classical  Adams-Bashforth-Moulton  method1^. 

Relative  error  buildup  is  controlled  by  an  error  tolerance  parameter,  so 
that  numerical  error  is  typically  kept  to  less  than  one  part  in  a  million 
for  times  up  to  10-50  picoseconds  during  a  simulation.  Symmetry  checks 
are  routinely  made  to  verify  the  accuracy  of  error  control. 

One  of  the  most  striking  features  of  the  CMD  runs  is  the  occurrence 
of  new  dimer  (and  also  higher-order)  molecular  arrangements  such  as 
[0:  :0]  as  can  be  clearly  seen  in  Figs.  4  and  6.  These  molecular 


arrangements  are  occurring  because  the  potentials  used  are  strictly 
central  and  of  the  additive  two-body  type  (i.e.,  the  form  of  an  OH 
potential  is  not  influenced  by  an  intervening  H  atom),  and  because  the  HH 
repulsion  is  weak  relative  to  the  strong  intermolecular  OH  interaction. 

The  new  water  arrangements  mentioned  in  the  above  paragraph  are  of 
course  not  realistic  in  that  they  are  ruled  out  by  the  molecular  orbital 
theory17  of  water-vapor  molecular  structure;  theoretically  the  OH 
potentials  should  not  be  the  simple  central  potentials  used  here,  but 
should  have  the  complicated  angular  dependence  characterizing  the 
lone-pair  electron  wavefunctions.17  Additionally,  the  Q  value  of 
0.139  e  resulting  from  Eq.  (3)  is  approximately  a  factor  of  1/3  the 
molecular  orbital  atomic-charge  population  prediction  of  Morokuma  and 
Pedersen11.  A  larger  Q  value  would  have  increased  the  intermolecular 
HH  repulsion  and  made  forming  such  new  molecular  arrangements  more 
difficult. 

At  100  kbar  (run  3,  Fig.  6)  the  shock  energy  associated  with  the  thin 

flyer  is  quickly  dissipated  in  breaking  the  intermolecular  OH  bonds  with 

the  result  that  a  volume  strain  does  not  propagate  into  the  sample  much 

beyond  21  ticks.  Yet,  because  of  the  energetically-favored 

above-mentioned  new  molecular  arrangements,  a  rotational  disturbance 

continues  to  propagate  to  the  end  of  the  target.  While  probably  not 

found  in  real  water,  such  rotational  propagation  effects  are  of 

considerable  interest  in  themselves  and  could  correspond  to  the  diffusion 

of  anisotropy  perturbations  (in  the  absence  of  strain)  as  treated  by 
18 

Leontovich  in  his  study  of  optical  scattering. 

At  the  300  kbar  level  of  runs  1  and  2  there  is  enough  strain  energy, 
even  with  the  thin  flyer,  for  the  volume  strain  to  propagate  to  the  end 
of  the  (short)  target  and  cause  spallation.  Additionally,  the  initial 
structures  of  runs  1  and  2  are  not  as  conducive  to  the  new  molecular 
arrangements  as  is  run  3.  The  result  is  that  in  runs  1  and  2  the  new 
arrangements  occur  only  after  sufficient  time  for  the  shock-induced 
temperature  to  result  in  molecular  orientation  relationships  favorable  to 
the  new  dimer  (and  higher-order)  phases. 

The  alternate  row  repeating  rotational  behavior  observed  in  runs  1 
and  2  is  the  result  of  impact  symmetry  and  the  already  mentioned  lack  of 
HH  interaction  across  the  rows  with  the  initial  triangular  structures. 
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While  in  run  1  this  repeating  behavior  has  largely  disappeared  by  the 
38  tick  frame,  it  is  still  very  much  in  evidence  in  run  2  at  36  ticks. 

The  reason  for  the  temporal  difference  between  runs  1  and  2  is  that  the 
Morse  HH  intermolecular  interactions  of  run  1  are  more  favorable  to  the 
new  molecular  arrangements  than  the  repulsive  Coulomb  interactions  of  run 
2.  The  result  is  that  these  new  molecular  arrangements  tend  to  form 
quickly  in  run  1,  and  in  the  process  of  forming  they  wash  out  the  every 
other  row  rotational  behavior. 

The  75  tick  run  3  frame  of  Fig.  6  is  very  interesting  in  that  it 
appears  that  the  high  degree  of  symmetry  perpendicular  to  the  direction 
of  shock  propagation  is  beginning  to  thermal ize  to  a  (probably)  more 
isotropic  state  characteristic  of  a  liquid.  This  effect  is  most  visible 
near  the  impact  line.  That  a  thermal ization  should  begin  to  occur  at 
67  ticks  (75  -  8)  is  in  agreement  with  our  present  understanding4  of  a 
picosecond  shock-front  rise  time  in  water  to  100  kbar.  While  this 
thermal ization  time  behavior  is  gratifying  it  should  be  viewed  with  some 
caution  because  our  CMD  structures  and  potentials  do  not  quite  duplicate 
real  water. 

The  intermolecular  00  potential  was  softened  (see  Table  III)  in  going 
from  run  1  to  run  2  in  order  to  lower  the  shock  velocity  to  a  value  more 

O 

in  keeping  with  experimental  observation  for  real  water.  While  it  is 
always  possible  for  a  given  structure  to  duplicate  bulk  shock-wave 
properties  by  varying  potential  parameters,  such  duplication  should  not 
be  taken  as  the  mark  of  success  in  that  many  materials  of  differing 
molecular  structure  have  very  similar  bulk  shock-wave  properties.  In 
that  regard  it  is  gratifying  to  observe  expected  changes  in  the 
intermolecular  HH  separation  in  all  the  runs.  However  for  the  formation 
of  complex  ionic  species  at  300  kbar  (runs  1  and  2),  it  is  possible  that 
we  may  require,  in  addition  to  a  pressure-dependent  lowering  of  the 
potential  energy  barrier  for  proton  transfer,  a  period  of  time  somewhat 
>  1  ps  in  order  to  have  sufficient  equipartitioning  of  energy  among  the 
molecular  degrees  of  freedom.  This  is  by  no  means  certain,  of  course, 
and  represents  an  area  for  further  study. 
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CONCLUSIONS 

Two-dimensional  hydrogen-bonded  structures  of  water-vapor  molecules 
have  been  impacted  under  CMO  conditions.  In  so  far  as  the  most  important 
aspect  of  liquid  water  (or  ice)  is  the  intermolecular  hydrogen  bonding, 
we  have  performed  CMD  shock-wave  experiments  on  water  (ice).  There 
still,  however,  remains  much  to  be  done.  Longer  running  times  are 
necessary  so  that  the  details  of  thermal ization  can  be  investigated  as  a 
function  of  time  and  so  that  a  supposedly  picosecond  shock-front 
rise-time  can  be  viewed  as  an  entity. 

In  order  to  have  a  model  that  can  yield  dissociation  at  300  kbar,  the 

0-H - 0  potentials  must  be  modified  by  first  allowing  for  the  presence 

of  a  double  minima,  then  secondly  permitting  under  shock  loading  the 
reduction  or  elimination  of  the  potential  barrier  inhibiting  proton 
transfer.  Thus  the  probability  of  forming  positive  ion  species  would  be 
enormously  enhanced.  If,  additionally,  the  simulations  are  then  run  to 
longer  times  in  order  to  ensure  thermal ization  (i.e.,  equipartition  of 
the  shock  energy  among  the  molecular  degrees  of  freedom)  and  complex 
formation  does  not  occur,  a  different  interpretation  of  the  saturating 
electrical  conductivity  may  be  necessary.  Finally,  the  potentials  must 
be  modified  to  be  more  representative  of  the  non-central  lone-pair 
electron  wavefunction  distribution.  That  latter  modification  will 
probably  be  sufficient  to  prevent  the  peculiar  molecular  arrangements 
discussed  in  the  previous  section.  While  the  various  PVT  phases  of  water 
are  probably  truly  3-dimensional  in  nature,  at  a  later  stage  something 
also  must  be  added  within  the  framework  of  our  2-dimensional  physics 
approach  to  allow  for  phases  which  correlate  with  those  of  real  water. 

The  small  intermolecular  HH  interaction  cut-off  distance  (2.05  A  in 
runs  1  and  2)  should  be  increased  to  some  value  which  is  understood  on 
the  basis  of  real  intermolecular  interactions.  If  Rc  were  thus 
increased  the  isolation  of  rotation  between  successive  molecular  rows 
would  be  diminished  with  perhaps  a  concurrent  faster  shock-front  rise 
time.  Such  a  lengthening  of  Rc  would  also  be  important  to  any  ultimate 
CMO  experiments  attempting  to  duplicate  shock-polarization  experiments. 

Much  work  remains  to  be  done. 
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Figure  1.  Initial  triangular  lattice  structure.  Large  circles  are 
oxygens  and  small  circles  are  hydrogens. 

Figure  2.  Initial  square  lattice  structure.  Large  circles  are  oxygens 
and  small  circles  are  hydrogens.  The  hydrogen  atom  labels 
Ha  and  Hb  are  for  purposes  of  discussion  within  the  text. 

Figure  3.  Coulomb  potential  for  runs  2  and  3.  The  cut-off  for  run  2  is 
Rc  =  2.05  A  and  for  run  3  is  Rc  =  2.3  A. 

Figure  4.  CMD  run  1.  Impact  occurs  at  approximately  8  ticks.  The 
impacting  plate  enters  from  the  left  and  is  composed  of  4 
columns  of  3  molecules  per  column.  Free  surface  spall  has 
already  occurred  at  26  ticks  with  the  spalled  layer  being  off 
screen  by  38  ticks.  Unconnected  hydrogens  are  not  unbonded, 
but  have  simply  been  reflected  across  the  transverse  periodic 
boundaries.  Notice  the  approximate  every  other  row  repeating 
pattern  of  rotational  behavior  to  26  ticks.  One  tick  = 

10"14  sec. 

Figure  5.  CMD  run  2.  Impact  occurs  at  approximately  7  ticks.  Again 

the  impacting  plate  enters  from  the  left  and  is  composed  of  4 
columns  of  3  molecules  per  column.  Spall  layer  separation  at 
the  free  surface  is  just  beginning  to  occur  at  36  ticks. 

Again  notice  the  approximate  every  other  row  repeating 
pattern  of  the  rotational  behavior,  but  now  to  36  ticks. 

One  tick  =  10-14  sec. 

Figure  6.  CMD  run  3.  Impact  is  just  occurring  at  8  ticks.  Impacting 
plate  enters  from  the  left  and  is  composed  of  4  columns  of  4 
molecules  per  column.  Notice  the  rotational  behavior  and  the 
peculiar  water  dimer  arrangements.  One  tick  *  10"^4  sec. 


Table  I 


SUMMARY  OF  WATER  IMPACT  CONFIGURATIONS 


RUN  NUMBER* 

STARTING 

LATTICE 

STRUCTURE 

IMPACT 

VELOCITY 
(cm/ sec) 

INTERMOLECULAR  POTENTIALS 

H-H  0-H  0-0 

1  (F3H20B) 

TRIANGULAR 

7.5  x  105 

MORSE 

MORSE 

MORSE 

2  (F3DF2C) 

TRIANGULAR 

7.5  x  105 

REPULSIVE 

COULOMB 

MORSE 

MORSE 

3  (F3SQUA) 

SQUARE 

4.0  x  105 

REPULSIVE 

COULOMB 

MORSE 

MORSE 
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Table  II 


INTRAMOLECULAR  MORSE  POTENTIAL  CONSTANTS* 


RUN  NO. 

INTERACTION 

De  (eV) 

A^A"1) 

Re  (A) 

Rc  (A) 

1 

OH 

5.1768 

2.2188 

1 .01 

3.5 

* 

HH 

0.5275 

2.1920 

1.01 

3.5 

2 

OH 

5.1768 

2.2188 

1.01 

3.5 

HH 

0.5275 

2.1920 

1.01 

3.5 

— -  r-* — — 

r  - 

3 

OH 

5.1768 

2.2188 

1.01 

2.3 

Vv 

HH 

0.5275 

2.1920 

1 .4284 

2.3 

*  Numbers  to  two  or  three  significant  figures  (such  as  2.3  and  1.01)  are  in 
reality  followed  by  zeros  (eg.,  2.30000...).  All  other  listed  numbers 
simply  represent  the  first  few  significant  figures. 

1  eV  =  1 .602  x  KT12  erg. 

a  Rc  is  the  cut-off  distance  (i.e.,  for  R  >  Rc  the  potential  is 
constant  and  the  force  vanishes). 


*  See  footnote  to  Table  II.  • 

a  See  corresponding  comment.  Table  II. 
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